/*! \file tvMosaicParam.cxx
 *      \brief Implements the methods to reconstruct 
 *       STP image imaged on TissueCyte 1000
*/
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include "tvMosaicParam.h"
#include <itkRawImageIO.h>
#include <itkImage.h>
//#include <itkCastImageFilter.h>
#include <itkRegionOfInterestImageFilter.h> 
#include <itkAdaptiveHistogramEqualizationImageFilter.h>
#include <itkFlipImageFilter.h>
#include <itkAddImageFilter.h>
#include <itkDivideImageFilter.h>
//#include <itkOrientImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkMultiplyByConstantImageFilter.h>
#include <itkPasteImageFilter.h>
#include <itkScaleTransform.h>
#include <itkIdentityTransform.h>
#include <itkResampleImageFilter.h>
//#include <omp.h>

typedef unsigned short TissueVisionPixelType;
typedef itk::Image<TissueVisionPixelType,2> TissueVisionImageType;
typedef itk::ImageFileReader<TissueVisionImageType> RawReaderType;

typedef unsigned short TilePixelType;
typedef itk::Image<TilePixelType,2> TileImageType;
typedef itk::Image<TilePixelType,2> AverageImageOutputType;
typedef itk::ImageFileReader<TileImageType> TileImageReaderType;
typedef itk::ImageFileWriter<TileImageType> TileWriterType;

typedef float AverageTilePixelType;
typedef itk::Image<AverageTilePixelType,2> AverageTileType;
typedef itk::ImageFileWriter<AverageImageOutputType> AverageImageWriterType;

typedef itk::CastImageFilter<TileImageType,AverageTileType> CastFilterType;
typedef itk::CastImageFilter<AverageTileType,AverageTileType> NullFilterType;
typedef itk::CastImageFilter<AverageTileType,TileImageType> AverageImageOutputFilterType;

typedef itk::DivideImageFilter <AverageTileType, AverageTileType, AverageTileType> DivideImageFilterType;

enum imageFormat {BIN_FORMAT = 0, TIF_FORMAT = 1};

using namespace std;
using namespace itk;

/*! \fn void tvMosaicParam::tvMosaicParam() 
 *  \brief Default Constructor
 *  \param None
 *  \exception None
 *  \return Reference to the current object
 */
tvMosaicParam::tvMosaicParam () {
    this->setDefaultParameters();
}

void tvMosaicParam::setDefaultParameters () {
  doConvertTiles = true;
  doCalculateCorrectionTile = true;
  doIlluminationCorrection = true;
  doStitching = true;
  doImageResize = true;

  doConvertTiles = false;
  doCalculateCorrectionTile = false;
  doIlluminationCorrection = false;
}

/*! \fn void tvMosaicParam::SetXOverlap(float xOverlap) 
 *  \brief Set the percentage of overlap along X-axis between tiles
 *  \param xOverlap Percentage of overlap between tiles (float)
 *  \exception None
 *  \return None
 */
void tvMosaicParam::SetXOverlap (float xOverlap) {
    xoverlap = xOverlap;
}

float tvMosaicParam::GetXOverlap () {
    return xoverlap;
}

float tvMosaicParam::GetYOverlap () {
    return yoverlap;
}

void tvMosaicParam::SetYOverlap (float yOverlap) {
    yoverlap = yOverlap;
}

void tvMosaicParam::SetZOverlap (float zOverlap) {
    zoverlap = zOverlap;
}

float tvMosaicParam::GetZOverlap () {
    return zoverlap;
}

void tvMosaicParam::SetConvertTiles(char* sConvertTiles) {
    string tmp = sConvertTiles;
    if (atoi(tmp.c_str()) == 1) {
        SetConvertTiles(true);
    }
    else {
        SetConvertTiles(false);
    }
}

void tvMosaicParam::SetCalculateCorrectionTile(char* sCalculateCorrectionTile) {
    string tmp = sCalculateCorrectionTile;
    if (atoi(tmp.c_str()) == 1) {
        SetCalculateCorrectionTile(true);
    }
    else {
        SetCalculateCorrectionTile(false);
    }
}

void tvMosaicParam::SetIlluminationCorrection(char* sIlluminationCorrection) {
    string tmp = sIlluminationCorrection;
    if (atoi(tmp.c_str()) == 1) {
        SetIlluminationCorrection(true);
    }
    else {
        SetIlluminationCorrection(false);
    }
}

void tvMosaicParam::SetStitching(char* sStitching) {
    string tmp = sStitching;
    if (atoi(tmp.c_str()) == 1) {
        SetStitching(true);
    }
    else {
        SetStitching(false);
    }
}

void tvMosaicParam::SetImageResize(char* sImageResize) {
    string tmp = sImageResize;
    if (atoi(tmp.c_str()) == 1) {
        SetImageResize(true);
    }
    else {
        SetImageResize(false);
    }
}

void tvMosaicParam::SetConvertTiles (bool bConvertTiles) {
    doConvertTiles = bConvertTiles;
}

void tvMosaicParam::SetCalculateCorrectionTile (bool bCalculateCorrectionTile) {
    doCalculateCorrectionTile = bCalculateCorrectionTile;
}

void tvMosaicParam::SetIlluminationCorrection (bool bIlluminationCorrection) {
    doIlluminationCorrection = bIlluminationCorrection;
}

void tvMosaicParam::SetStitching (bool bStitching) {
    doStitching = bStitching;
}

void tvMosaicParam::SetImageResize (bool bImageResize) {
    doImageResize = bImageResize;
}

void tvMosaicParam::SetParamFileName (char* paramFileName) {
  tvParamFileName = paramFileName;
}

void tvMosaicParam::SetDataHome (char* dataSetHomeDirectory) {
  tvDataHome = dataSetHomeDirectory;
}

void tvMosaicParam::SetDataProcessingHome (char *dataSetProcessingHome) {
  tvDataProcessingHome = dataSetProcessingHome;
}

void tvMosaicParam::SetDataPrefix (char* dataSetPrefix) {
  tvDataPrefix = dataSetPrefix;
}

void tvMosaicParam::SetOutputImageFormat (char* outputImageFormat) {
  tvOutputImageFormat = outputImageFormat;
}

void tvMosaicParam::GenerateInputDirectory () {
  tvInputDirectory = tvDataHome + tvDataPrefix + "bin";
}

void tvMosaicParam::GenerateOutputDirectory() {
  tvOutputDirectory = tvDataProcessingHome + tvDataPrefix + tvOutputImageFormat;
}

void tvMosaicParam::GenerateCorrectedTileDirectory() {
  tvCorrectedTileOutputDirectory = tvDataProcessingHome + tvDataPrefix + "corrected_" + tvOutputImageFormat;
}

void tvMosaicParam::GenerateStitchedDirectory() {
  tvStitchedImageDirectory = tvDataProcessingHome + tvDataPrefix + "stitched/";
}

void tvMosaicParam::GenerateWarpingDirectory() {
  tvWarpedImageDirectory = tvDataProcessingHome + tvDataPrefix + "warping/";
}

void tvMosaicParam::GenerateCorrectionParameters () {
  tvAverageImageDirectory = tvDataProcessingHome + "Correction/";
  tvAverageImageFileName = tvAverageImageDirectory + tvDataPrefix + "AverageImage.tif";
}

void tvMosaicParam::runReconstructionPipeline() {
    this->initializeDataset();
    this->readMosaicParamFile();
    if (doConvertTiles) {
        this->writeImageTiles();
    }
    if (doCalculateCorrectionTile) {
        this->createAverageTile();
    }
    if (doIlluminationCorrection) {
        this->writeCorrectedImageTiles();
    }
    if (doStitching) {
        this->writeAllStitchedImages();
    }
    if (doImageResize) {
        this->writeScaledImage(0.05);
    }
}

void tvMosaicParam::initializeExampleDataset() {
  this->SetParamFileName("/Users/kannanuv/Documents/workspace/data/120120_babyset_ch3/Mosaic1_20_2012_4_20_35_PM_ch3.txt");
  this->SetDataHome("/Users/kannanuv/Documents/workspace/data/120120_babyset_ch3/");
  this->SetDataProcessingHome("/Users/kannanuv/Documents/workspace/data/120120_babyset_ch3/");
  this->SetDataPrefix("120120_babyset_ch3_");
  this->SetOutputImageFormat("tif");
  ifstream paramFile(tvParamFileName.c_str());
  if (!paramFile) {
    cout << "Cannot open example dataset parameterfile/ it doesn't exist:" << tvParamFileName << endl;
    exit (1);
  }
}

void tvMosaicParam::initializeDataset() {
  this->GenerateInputDirectory();
  this->GenerateOutputDirectory();
  this->GenerateCorrectedTileDirectory();
  this->GenerateStitchedDirectory();
  this->GenerateWarpingDirectory();

  croppedTileSizeX = 812;
  croppedTileSizeY = 812;
  croppedStartX = 10;
  croppedStartY = 15;

  this->SetXOverlap(0.12);
  this->SetYOverlap(0.12);
  this->GenerateCorrectionParameters();
}

int tvMosaicParam::readMosaicParamFile () {
  cout << "Reading parameter file... ";

  // open a text file for reading
  ifstream paramFile(tvParamFileName.c_str());
  char str[255];

  if (!paramFile) {
    cout << "Cannot open file:" << tvParamFileName << endl;
	exit (1);
  }

  printf("File opened. Now parsing it...\n");
  while (paramFile) {
    paramFile.getline (str, 255);
    //cout << str << endl;
    string line(str);
	string keyName;
	string keyValue;
	size_t tokenPos;
	size_t endOfLine;
	tokenPos = line.find_first_of (":");
	endOfLine = line.find_first_of ("\n");
	keyName.assign (str, tokenPos);
	keyValue.assign (str, tokenPos+1, endOfLine-tokenPos);
    //cout << "Token Position = " << tokenPos << "Key Name : " << keyName << "Key value : " << keyValue << endl;

    if (keyName.compare ("rows") == 0) {
	  rows = atoi (keyValue.c_str());
	} else if (keyName.compare ("columns") == 0) {
	  columns = atoi (keyValue.c_str());
	} else if (keyName.compare ("layers") == 0) {
	  layers = atoi (keyValue.c_str());
	} else if (keyName.compare ("xres") == 0) {
	  xres = atoi (keyValue.c_str());
	} else if (keyName.compare ("yres") == 0) {
	  yres = atoi (keyValue.c_str());
	} else if (keyName.compare ("zres") == 0) {
	  zres = atoi (keyValue.c_str());
	} else if (keyName.compare ("mrows") == 0) {
	  mrows = atoi (keyValue.c_str());
	} else if (keyName.compare ("mcolumns") == 0) {
	  mcolumns = atoi (keyValue.c_str());
	} else if (keyName.compare ("mrowres") == 0) {
	  mrowres = atoi (keyValue.c_str());
	} else if (keyName.compare ("mcolumnres") == 0) {
	  mcolumnres = atoi (keyValue.c_str());
	} else if (keyName.compare ("sections") == 0) {
	  sections = atoi (keyValue.c_str());
	} else if (keyName.compare ("sectionres") == 0) {
	  sectionres = atoi (keyValue.c_str());
	} else if (keyName.compare ("channels") == 0) {
	  channels = atoi (keyValue.c_str());
	} else if (keyName.compare ("Pixrestime") == 0) {
	  Pixrestime = atoi (keyValue.c_str());
	} else if (keyName.compare ("xoverlap") == 0) {
	  xoverlap = atoi (keyValue.c_str());
	} else if (keyName.compare ("yoverlap") == 0) {
	  yoverlap = atoi (keyValue.c_str());
	} else if (keyName.compare ("zoverlap") == 0) {
	  zoverlap = atoi (keyValue.c_str());
	} else if (keyName.compare ("normflag") == 0) {
	  normflag = atoi (keyValue.c_str());
	} else if (keyName.compare ("Normimage") == 0) {
	  Normimage = atoi (keyValue.c_str());
	} else if (keyName.compare ("startnum") == 0) {
	  startnum = atoi (keyValue.c_str());
	} else if (keyName.compare ("excwavelength") == 0) {
	  excwavelength = atoi (keyValue.c_str());
	} else if (keyName.compare ("excpower") == 0) {
	  excpower = atoi (keyValue.c_str());
	} else if (keyName.compare ("numfoci") == 0) {
	  numfoci = atoi (keyValue.c_str());
	// String attributes TBD
	} else if (keyName.compare ("objective") == 0) {
	  objective = keyValue;
	} else if (keyName.compare ("basefilename") == 0) {
	  basefilename = keyValue;
	} else if (keyName.compare ("fSuffix") == 0) {
	  fSuffix = keyValue;
	} else if (keyName.compare ("fPrefix") == 0) {
	  fPrefix = keyValue;
	} else if (keyName.compare ("sample") == 0) {
	  sample = keyValue;
	} else if (keyName.compare ("acqDate") == 0) {
	  acqDate = keyValue;
	}

    // Comments is the last of the sting attributes then starts
	// the postion information of each tile
    if (keyName.compare ("comments") == 0) {
	  cout << "Reading parameters... Done" << endl;
	  break;
	}
  }
  // Get info about the tile positions
  string displacement;
  string position;
  stringstream ss;
  stringstream ss1;
  size_t tokenPos;
  size_t endOfLine;
  int iTileX = 1;
  int iTileY = 1;
  int iTileZ = -1;
  int tileZID = 0;
  int displacementX = 0;
  int displacementY = 0;
  string fileName1;
  string fileName2;
  string fPrefix1("");
  string fSuffix1("");
  fPrefix1.append(fPrefix,0,fPrefix.length()-1);
  fSuffix1.append(fSuffix,0,fSuffix.length()-1);
  //char *cDataOutputFormat;
  //cDataOutputFormat = tvOutputImageFormat.c_str();
  for (int iTile = 0; iTile < mrows * mcolumns * sections; iTile++) {
	// X Position
    paramFile.getline (str, 255);
    string Xline(str);
    tokenPos = Xline.find_first_of (":");
    endOfLine = Xline.find_first_of ("\n");
    displacement.assign (str, 4, tokenPos-4);
    position.assign (str, tokenPos+1, endOfLine-tokenPos);
    XPos.push_back (atoi (position.c_str()));
    displacementX = atoi (displacement.c_str());
    if (displacementX > 0) {
      iTileX = iTileX + 1;
    } else if (displacementX < 0){
      iTileX = iTileX - 1;
    }
	  
    // Y Position
    paramFile.getline (str, 255);
    //cout << "Line " << ": " << str << endl;
    string Yline(str);
    tokenPos = Yline.find_first_of (":");
    endOfLine = Yline.find_first_of ("\n");
    displacement.assign (str, 4, tokenPos-4);
    position.assign (str, tokenPos+1, endOfLine-tokenPos);
    YPos.push_back (atoi (position.c_str()));
    displacementY = atoi (displacement.c_str());
    if (displacementY > 0) {
      iTileY = iTileY + 1;
    } else if (displacementY < 0){
      iTileY = iTileY - 1;
    }
    displacement.clear();
    position.clear();

    if (displacementX == 0 && displacementY == 0) {
       iTileZ = iTileZ + 1;
       iTileY = 1;
    }
    for (int iLayer = 0; iLayer < layers; iLayer++) {
      tileX.push_back (iTileY);
      tileY.push_back (iTileX);
      tileZID = iTileZ * layers + iLayer;
      tileZ.push_back (tileZID);
  
      // I don't know why I am doing this. It doesn't work unless I use
      // prefix1 & fSuffix. May be trailing white spaces.
      ss1.str("");
      //ss1 << fPrefix1 << " " << startnum + (iTile*layers)+iLayer << " " << fSuffix1;
      
      ss1 << fPrefix1 << " " << startnum + (iTile*layers)+iLayer << fSuffix1;
	  fileName1 = ss1.str();

      ss << "Tile_Z" << setfill ('0') << setw(3) << tileZID 
         << "_Y" << setw(3) << iTileX
         << "_X" << setw(3) << iTileY << "." << tvOutputImageFormat.c_str();// << endl;
      fileName2 = ss.str();
      ss.str("");
      binaryFileName.push_back (fileName1);
      imageFileName.push_back (fileName2);
      layerID.push_back (iLayer);
    }
  }
  paramFile.close ();

  return 0;
}

//*****************************************************************************
// int tvMosaicParam::writeImageTiles () 
//*****************************************************************************
// converts all the binary files to image files
//*****************************************************************************
int tvMosaicParam::writeImageTiles () {
  unsigned int dimension = 2;
  string dataOutputFormat(tvOutputImageFormat.c_str());
  string rawImageDir(tvInputDirectory);
  string outputImageDir(tvOutputDirectory);
  typedef itk::CastImageFilter<TissueVisionImageType, TileImageType> CastFilterType; 
  typedef itk::RegionOfInterestImageFilter<TissueVisionImageType, TissueVisionImageType> CropFilterType;
  typedef itk::FlipImageFilter<TileImageType> ImageFlipperType;
  ofstream stdOutLog;
  ofstream stdErrLog;
  stdOutLog.open ("StdOut.log");
  stdErrLog.open ("StdErr.log");

  // Wrting log files
  stdOutLog << "Raw Images directory: " << rawImageDir << endl;
  stdOutLog << "Output image directory: " << outputImageDir << endl;

  TissueVisionImageType::SizeType inputSize;
  inputSize[0] = rows;
  inputSize[1] = columns;
  cout << "Rows = " << rows << ", Columns = " << columns << endl;
  TissueVisionImageType::SizeType outputSize;
  outputSize[0] = croppedTileSizeX;
  outputSize[1] = croppedTileSizeY;
  cout << "Cropped Rows = " << croppedTileSizeX << ", Cropped Columns = " << croppedTileSizeY << endl;
  int imageFormat = TIF_FORMAT;
  RawReaderType::Pointer reader = RawReaderType::New();

  //omp_set_num_threads(8);
  //#pragma omp parallel for //UVK TBD
  //for (int iTile = 0; iTile < 10; iTile++) {
  for (int iTile = 0; iTile < (int) binaryFileName.size(); iTile++) {
    string *binFileName = new string("");
    string *imFileName = new string("");

    if (imageFormat == BIN_FORMAT) {
      // Set RAW image reader properties
      itk::RawImageIO<TilePixelType, 2>::Pointer rawImageIO
       = itk::RawImageIO<TilePixelType, 2>::New();
      rawImageIO->SetByteOrderToLittleEndian();
      rawImageIO->SetFileDimensionality (dimension);
      rawImageIO->SetDimensions (0, rows);
      rawImageIO->SetDimensions (1, columns);
      rawImageIO->SetNumberOfComponents(1); 
      reader->SetImageIO (rawImageIO);
    } 
    else {
      //Do nothing
    }

    // Crop filter instantiations
    TissueVisionImageType::RegionType inputRegion;
    TissueVisionImageType::RegionType::IndexType inputStart;
    inputStart[0] = croppedStartX;
    inputStart[1] = croppedStartY;
    inputRegion.SetIndex(inputStart);
    const TissueVisionImageType::SpacingType& spacing = reader->GetOutput()->GetSpacing();
    const TissueVisionImageType::PointType& inputOrigin = reader->GetOutput()->GetOrigin();
    double  outputOrigin [2];
    inputRegion.SetSize (outputSize);

    for(unsigned int i=0; i<2 ; i++) {
      outputOrigin[i] = inputOrigin[i] + spacing[i] * inputStart[i];
    }
  
    CropFilterType::Pointer cropFilter = CropFilterType::New();
    cropFilter->SetRegionOfInterest(inputRegion);
  
    CastFilterType::Pointer caster = CastFilterType::New(); 
  
    TileWriterType::Pointer writer = TileWriterType::New();
  
    binFileName->clear();
    binFileName->append(rawImageDir);
    binFileName->append("/");
    binFileName->append(binaryFileName.at(iTile).c_str());
    imFileName->clear();
    imFileName->append(outputImageDir);
    imFileName->append("/");
    imFileName->append(imageFileName.at(iTile).c_str());
 
    //cout << "Binary File: " << binFileName->c_str() << endl;
    reader->SetFileName (binFileName->c_str());

    // Cropping to the edge pixels
    cropFilter->SetInput (reader->GetOutput());
    cropFilter->GetOutput()->SetSpacing (spacing);
    cropFilter->GetOutput()->SetOrigin (outputOrigin);
 
    // Casting data from input pixel type to output pixel type
    caster->SetInput (cropFilter->GetOutput());

    // Flip Image Filter
    ImageFlipperType::Pointer imageFlipper = ImageFlipperType::New();
    ImageFlipperType::FlipAxesArrayType flipAxes;
    if (imageFormat == BIN_FORMAT){
      flipAxes[0] = false;
      flipAxes[1] = true;
    }
    else {
      flipAxes[0] = false;
      flipAxes[1] = false;
    }
    imageFlipper->SetFlipAxes(flipAxes);
    imageFlipper->SetInput (caster->GetOutput());

    //Writer setup
    writer->SetFileName (imFileName->c_str());
    writer->SetInput (imageFlipper->GetOutput());
    try {
      writer->Update(); 
      stdOutLog << "Binary File: " << binFileName->c_str() << endl;
      stdOutLog << "Tile ID: " << iTile << ", Image File: " << imFileName->c_str() << endl;
      //cout << "Tile ID: " << iTile << ", Image File: " << imFileName->c_str() << endl;
    }
    catch (itk::ExceptionObject & e) {
      std::cerr << "Exception in processing file" << std::endl;
      std::cerr << "Binary File: " << binFileName->c_str() << endl;
      std::cerr << "Tile ID: " << iTile << ", Image File: " << imFileName->c_str() << endl;
      std::cerr << e << std::endl;

      stdErrLog << "Exception in processing file" << std::endl;
      stdErrLog << "Binary File: " << binFileName->c_str() << endl;
      stdErrLog << "Tile ID: " << iTile << ", Image File: " << imFileName->c_str() << endl;
      stdErrLog << e << std::endl;
    }
  }

  stdOutLog.close();
  stdErrLog.close();
  return 0;
}

void tvMosaicParam::CreateConstantImage(TileImageType::Pointer image, TilePixelType pixelValue)
{
  // Create an image with 2 connected components
  TileImageType::IndexType start;
  start[0] = 0;
  start[1] = 0;

  TileImageType::SizeType size;
  size[0] = croppedTileSizeX;
  size[1] = croppedTileSizeY;

  TileImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);

  image->SetRegions(region);
  image->Allocate();
  image->FillBuffer(pixelValue);
}

/*! \fn void tvMosaicParam::CreateConstantImage(AverageTileType::Pointer image, AverageTilePixelType pixelValue)
 *  \brief Utility function to create an constant image
 *  \param image Image pointer
 *  \param pixelValue Intensity value set to all pixels of the image
 *  \exception None
 *  \return None
 */
void tvMosaicParam::CreateConstantImage(AverageTileType::Pointer image, AverageTilePixelType pixelValue)
{
  // Create an image with 2 connected components
  TileImageType::IndexType start;
  start[0] = 0;
  start[1] = 0;

  TileImageType::SizeType size;
  size[0] = croppedTileSizeX;
  size[1] = croppedTileSizeY;

  TileImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);

  image->SetRegions(region);
  image->Allocate();

  image->FillBuffer(pixelValue);
}

int tvMosaicParam::createAverageTile() {
    cout << "Creating correction tile... ";
  typedef itk::AddImageFilter<AverageTileType, AverageTileType, AverageTileType> AddImageFilterType;

  TileImageType::SizeType outputSize;
  outputSize[0] = croppedTileSizeX;
  outputSize[1] = croppedTileSizeY;

  AddImageFilterType::Pointer addImageFilter = AddImageFilterType::New();
  DivideImageFilterType::Pointer divideImageFilter = DivideImageFilterType::New();
  AverageImageWriterType::Pointer averageImageWriter = AverageImageWriterType::New();

  // Setting up blank image as average
  NullFilterType::Pointer nullFilter = NullFilterType::New();
  TileImageReaderType::Pointer reader = TileImageReaderType::New();
  CastFilterType::Pointer castFilter = CastFilterType::New(); 
  AddImageFilterType::Pointer addFilter = AddImageFilterType::New();
  AverageImageOutputFilterType::Pointer averageImageOutputFilterType = AverageImageOutputFilterType::New();
  AverageTileType::Pointer divisorImage = AverageTileType::New();
  AverageTileType::Pointer dividendImage = AverageTileType::New();

  CreateConstantImage (divisorImage, (float)imageFileName.size());
  CreateConstantImage (dividendImage, (float)0.0);

  // Set up the itk pipeline
  string *imFileName = new string("");
  string outputImageDir(tvOutputDirectory);
  for (int iTile = 0; iTile < (int) imageFileName.size(); iTile++) {
  //for (int iTile = 0; iTile < 10; iTile++) {
    //cout << "Processing iTile = " << iTile << ", " << imFileName->c_str() << endl;
    imFileName->clear();
    imFileName->append(outputImageDir);
    imFileName->append("/");
    imFileName->append(imageFileName.at(iTile).c_str());
    reader->SetFileName (imFileName->c_str());
    castFilter->SetInput (reader->GetOutput());
    castFilter->Update();
    addFilter->SetInput1 (castFilter->GetOutput());
    addFilter->SetInput2 (dividendImage);
    dividendImage->SetOrigin (castFilter->GetOutput()->GetOrigin());
    dividendImage->SetSpacing (castFilter->GetOutput()->GetSpacing());
    nullFilter->SetInput (addFilter->GetOutput());
    dividendImage = nullFilter->GetOutput();
    nullFilter->Update();
  }

  divideImageFilter->SetInput1(dividendImage);
  divideImageFilter->SetInput2(divisorImage);
  divisorImage->SetOrigin (dividendImage->GetOrigin());
  divisorImage->SetSpacing (dividendImage->GetSpacing());
  averageImageOutputFilterType->SetInput (divideImageFilter->GetOutput());
  averageImageWriter->SetInput(averageImageOutputFilterType->GetOutput());

  averageImageWriter->SetFileName (tvAverageImageFileName.c_str());
  cout << "Done" << endl << "Average Image name : " <<  tvAverageImageFileName.c_str() << endl;
  averageImageWriter->Update();
  return 0;
}

int tvMosaicParam::writeCorrectedImageTiles() {
  cout << "tvMosaicParam::writeCorrectedImageTiles" << imageFileName.size() << endl;

  typedef MinimumMaximumImageCalculator <TileImageType> ImageCalculatorType;
  ImageCalculatorType::Pointer imageCalculator = ImageCalculatorType::New();

  typedef MultiplyByConstantImageFilter <AverageTileType,float,AverageTileType> MultiplyByConstantImageType;
  MultiplyByConstantImageType::Pointer imageMultiplier = MultiplyByConstantImageType::New();

  string *inputTileName = new string("");
  string *correctedTileName = new string("");
  string inputImageDir(tvOutputDirectory);
  string outputImageDir(tvCorrectedTileOutputDirectory);
  TileImageReaderType::Pointer tileReader = TileImageReaderType::New();
  CastFilterType::Pointer tileCastFilter = CastFilterType::New(); 
  AverageImageOutputFilterType::Pointer averageImageOutputFilterType = AverageImageOutputFilterType::New();
  TileWriterType::Pointer correctedTileWriter = TileWriterType::New();

  //Read the average tile
  TileImageReaderType::Pointer averageTileReader = TileImageReaderType::New();
  AverageTileType::Pointer averageTile = AverageTileType::New();
  CastFilterType::Pointer castFilter = CastFilterType::New(); 
  CreateConstantImage (averageTile, 0.0);
  averageTileReader->SetFileName (tvAverageImageFileName.c_str());
  averageTileReader->Update();
  castFilter->SetInput (averageTileReader->GetOutput());
  averageTile = castFilter->GetOutput();

  DivideImageFilterType::Pointer divideImageFilter = DivideImageFilterType::New();

  //Get average tile image minimum
  imageCalculator->SetImage(averageTileReader->GetOutput());
  imageCalculator->ComputeMinimum();
  imageMultiplier->SetConstant((float) imageCalculator->GetMinimum());

  for (int iTile = 0; iTile < (int) imageFileName.size(); iTile++) {
    //cout << "Processing iTile = " << iTile;
    inputTileName->clear();
    inputTileName->append(inputImageDir);
    inputTileName->append("/");
    inputTileName->append(imageFileName.at(iTile).c_str());
    //cout << ", Image name : " <<  inputTileName->c_str() << endl;
    tileReader->SetFileName (inputTileName->c_str());
    tileCastFilter->SetInput (tileReader->GetOutput());
    divideImageFilter->SetInput1(tileCastFilter->GetOutput());
    divideImageFilter->SetInput2(averageTile);
    averageTile->SetOrigin (tileCastFilter->GetOutput()->GetOrigin());
    averageTile->SetSpacing (tileCastFilter->GetOutput()->GetSpacing());


    imageMultiplier->SetInput(divideImageFilter->GetOutput());
    imageMultiplier->Update();

    averageImageOutputFilterType->SetInput(imageMultiplier->GetOutput());
    correctedTileWriter->SetInput(averageImageOutputFilterType->GetOutput());
    correctedTileName->clear();
    correctedTileName->append(outputImageDir);
    correctedTileName->append("/");
    correctedTileName->append(imageFileName.at(iTile).c_str());
    correctedTileWriter->SetFileName(correctedTileName->c_str());
    correctedTileWriter->Update();
  }
  return 0;
}

string tvMosaicParam::getStitchedImageFileName (unsigned int sectionID, unsigned int layerID) {
    stringstream mosaicFileNameStream;
    stringstream mosaicFullFileNameStream;
    mosaicFileNameStream.str();
    mosaicFullFileNameStream.str();

    mosaicFileNameStream << "StitchedImage"
                         << setfill('0') << setw(4) << sectionID;
    if (layers > 1) {
      mosaicFileNameStream << "_layer" << setw(3) << layerID << "." << tvOutputImageFormat.c_str();
    }
    stitchedFileName.push_back(mosaicFileNameStream.str());

    mosaicFullFileNameStream << tvStitchedImageDirectory << mosaicFileNameStream.str();
    return mosaicFullFileNameStream.str();
}


void tvMosaicParam::writeStitchedImageITK(unsigned int sectionID, unsigned int layerID) {
  unsigned int startTileIndex = 0;
  unsigned int endTileIndex = 0;
  unsigned int tilesPerSection = 0;
  tilesPerSection = mrows * mcolumns * layers;
  startTileIndex = sectionID * tilesPerSection;
  endTileIndex = startTileIndex + tilesPerSection - 1;

  // Create blank mosaic image
  TileImageType::IndexType mosaicIndex;
  TileImageType::IndexType mosaicIndex1;
  mosaicIndex[0] = 0;
  mosaicIndex[1] = 0;
  TileImageType::SizeType mosaicSize;
  mosaicSize[0] = mcolumns * croppedTileSizeX;
  mosaicSize[1] = mrows * croppedTileSizeY;
  TileImageType::RegionType region;
  region.SetSize(mosaicSize);
  region.SetIndex(mosaicIndex);

  TileImageType::Pointer mosaicImage = TileImageType::New();
  mosaicImage->SetRegions(region);
  mosaicImage->Allocate();
  mosaicImage->FillBuffer(40);

  //typedef itk::TranslationTransform<double,2> TranslationTransformType;
  //TranslationTransformType::Pointer transform = TranslationTransformType::New();
  //TranslationTransformType::OutputVectorType translation;
  //typedef itk::ResampleImageFilter<ImageType, ImageType> ResampleImageFilterType;
  //ResampleImageFilterType::Pointer resampleFilter = ResampleImageFilterType::New();

  TileImageReaderType::Pointer reader = TileImageReaderType::New();
  TileImageType::Pointer tileImage = TileImageType::New();
  typedef itk::PasteImageFilter <TileImageType, TileImageType> PasteImageFilterType;
  PasteImageFilterType::Pointer pasteFilter = PasteImageFilterType::New ();

  for (unsigned int iTile = startTileIndex + layerID; iTile < endTileIndex; iTile = iTile + layers) {
      string *imFileName = new string("");
      string outputImageDir(tvOutputDirectory);
      imFileName->clear();
      imFileName->append(outputImageDir);
      imFileName->append("/");
      imFileName->append(imageFileName.at(iTile).c_str());
      reader->SetFileName (imFileName->c_str());
      tileImage = reader->GetOutput();
      reader->Update();
      pasteFilter->SetSourceImage(tileImage);
      pasteFilter->SetDestinationImage (mosaicImage);
      pasteFilter->SetSourceRegion(tileImage->GetLargestPossibleRegion());

      mosaicIndex1[0] = (tileX.at(iTile) - 1) * croppedTileSizeX;
      mosaicIndex1[1] = (tileY.at(iTile) - 1) * croppedTileSizeY;
      pasteFilter->SetDestinationIndex(mosaicIndex1);
      mosaicImage = pasteFilter->GetOutput();
      region.SetIndex(mosaicIndex);
      mosaicImage->SetRegions(region);
      //translation[0] = -iX * 812;
      //translation[1] = -iY * 812;
      //transform->Translate(translation);
      //resampleFilter->SetTransform(transform.GetPointer());
      pasteFilter->Update();
  }

  cout << "Writing Mosaic Image" << endl;
  TileWriterType::Pointer mosaicImageWriter = TileWriterType::New();
  mosaicImageWriter->SetFileName(getStitchedImageFileName(sectionID, layerID).c_str());
  mosaicImageWriter->SetInput(mosaicImage);
  mosaicImageWriter->Update();  
}

void tvMosaicParam::writeAllStitchedImages() {
  for (unsigned int iSection = 0; iSection < sections; iSection++) {
    for (unsigned int iLayer = 0; iLayer < layers; iLayer++) {
      writeStitchedImage(iSection, iLayer);
    }
  }
}

void tvMosaicParam::writeStitchedImage(unsigned int sectionID, unsigned int layerID) {
  unsigned int startTileIndex = 0;
  unsigned int endTileIndex = 0;
  unsigned int tileUpperLeftX = 0;
  unsigned int tileUpperLeftY = 0;
  unsigned int tileLowerRightX = 0;
  unsigned int tileLowerRightY = 0;
  unsigned int tilesPerSection = 0;
  tilesPerSection = mrows * mcolumns * layers;
  startTileIndex = sectionID * tilesPerSection;
  endTileIndex = startTileIndex + tilesPerSection - 1;
  tileLowerRightX = croppedTileSizeX;
  tileLowerRightY = croppedTileSizeY;

  // Create blank mosaic image
  TileImageType::IndexType mosaicIndex;
  TileImageType::IndexType tileIndex;
  mosaicIndex[0] = 0;
  mosaicIndex[1] = 0;
  TileImageType::SizeType mosaicSize;
  mosaicSize[0] = mcolumns * croppedTileSizeX;
  mosaicSize[1] = mrows * croppedTileSizeY;
  TileImageType::RegionType region;
  region.SetSize(mosaicSize);
  region.SetIndex(mosaicIndex);

  TileImageType::Pointer mosaicImage = TileImageType::New();
  mosaicImage->SetRegions(region);
  mosaicImage->Allocate();
  mosaicImage->FillBuffer(0);

  TileImageReaderType::Pointer reader = TileImageReaderType::New();
  TileImageType::Pointer tile = TileImageType::New();

  for (unsigned int iTile = startTileIndex + layerID; iTile < endTileIndex; iTile = iTile + layers) {
      // Read the tile
      string *imFileName = new string("");
      string outputImageDir(tvOutputDirectory);
      imFileName->clear();
      imFileName->append(outputImageDir);
      imFileName->append("/");
      imFileName->append(imageFileName.at(iTile).c_str());
      reader->SetFileName (imFileName->c_str());
      //cout << imFileName->c_str() << endl;
      tile = reader->GetOutput();
      reader->Update();

      for (unsigned int iTileX = tileUpperLeftX; iTileX < tileLowerRightX; iTileX++) {
          for (unsigned int iTileY = tileUpperLeftY; iTileY < tileLowerRightY; iTileY++) {
              // Set tile and mosaic pixel co-ordinates
              tileIndex[0] = iTileX;
              tileIndex[1] = iTileY;
              mosaicIndex[0] = (tileX.at(iTile) - 1) * croppedTileSizeX + tileIndex[0];
              mosaicIndex[1] = (tileY.at(iTile) - 1) * croppedTileSizeY + tileIndex[1];
              mosaicImage->SetPixel (mosaicIndex, tile->GetPixel(tileIndex));
          }
      }
  }

  stringstream mosaicFileNameStream;
  stringstream mosaicFullFileNameStream;
  mosaicFileNameStream.str();
  mosaicFullFileNameStream.str();

  mosaicFileNameStream << "StitchedImage"
                       << setfill('0') << setw(4) << sectionID;
  if (layers > 1) {
    mosaicFileNameStream << "_layer" << setw(3) << layerID << "." << tvOutputImageFormat.c_str();
  }
  stitchedFileName.push_back(mosaicFileNameStream.str());

  mosaicFullFileNameStream << tvStitchedImageDirectory << mosaicFileNameStream.str();

  TileWriterType::Pointer mosaicImageWriter = TileWriterType::New();

  mosaicImageWriter->SetFileName(mosaicFullFileNameStream.str().c_str());
  cout << "Writing Mosaic Image" << mosaicFullFileNameStream.str().c_str() << endl;
  mosaicImageWriter->SetInput(mosaicImage);
  mosaicImageWriter->Update();

  /*for (unsigned int iColumn = 0; iColumn < mosaicSize[0]; iColumn++) {
      for (unsigned int iRow = 0; iRow < mosaicSize[1]; iRow++) {
          mosaicIndex[0] = iRow;
          mosiacIndex[1] = iColumn;
          for (unsigned int iTile = startTileIndex; iTile < endTileIndex; iTile++) {
              // check for contributing tiles to the pixel
              // blend the values of the pixels
              // write the pixel values to mosaic image
          }
      }
  }*/
}

void tvMosaicParam::writeScaledImage(float scaleFactor) {
  typedef itk::ResampleImageFilter<TileImageType, TileImageType> ResampleImageFilterType;
  typedef itk::IdentityTransform<double, 2> TransformType;
  TransformType::Pointer identityTransform = TransformType::New();
  TileImageReaderType::Pointer reader = TileImageReaderType::New();
  TileWriterType::Pointer writer = TileWriterType::New();
  TileImageType::SizeType inputSize;
  TileImageType::SizeType outputSize;
  TileImageType::SpacingType outputSpacing;

  ResampleImageFilterType::Pointer resampleFilter = ResampleImageFilterType::New();

  for(vector<string>::iterator it = stitchedFileName.begin(); it != stitchedFileName.end(); ++it) {
    string stitchedFileName = *it;
    string scaledImageName = tvWarpedImageDirectory + stitchedFileName;
    stitchedFileName = tvStitchedImageDirectory + stitchedFileName;
    reader->SetFileName(stitchedFileName.c_str());
    reader->Update();
    inputSize = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
    outputSize[0] = ceil (inputSize[0] * scaleFactor);
    outputSize[1] = ceil (inputSize[1] * scaleFactor);
    outputSpacing[0] = reader->GetOutput()->GetSpacing()[0] * (1 / scaleFactor);
    outputSpacing[1] = reader->GetOutput()->GetSpacing()[1] * (1 / scaleFactor);
    resampleFilter->SetInput(reader->GetOutput());
    resampleFilter->SetOutputSpacing(outputSpacing);
    resampleFilter->SetSize(outputSize);
    resampleFilter->SetTransform(identityTransform);
    cout << "Original image : " << stitchedFileName << endl;
    cout << "Down sizing image : " << scaledImageName << endl;

    writer->SetFileName(scaledImageName.c_str());
    writer->SetInput(resampleFilter->GetOutput());
    writer->Update();
  }
}
